Correlations and invariance of seismicity under renormalization- group transformations 



Alvaro Corral 

Grup de Ffsica Estadistica, Departament de Fi'sica, Facultat de Ciencies, 
Universitat Autdnoma de Barcelona, E-08193 Bellaterra, Barcelona, Spaing 
(Dated: February 2, 2008) 

The effect of transformations analogous to those of the real-space renormalization group are analyzed for 
the temporal occurrence of earthquakes. The distribution of recurrence times turns out to be invariant under 
such transformations, for which the role of the correlations between the magnitudes and the recurrence times 
are fundamental. A general form for the distribution is derived imposing only the self-similarity of the process, 
which also yields a scaling relation between the Gutenberg-Richter b— value, the exponent characterizing the 
correlations, and the recurrence-time exponent. This approach puts the study of the structure of seismicity in 
the context of critical phenomena. 



The study of statistical seismology has a long history, ex- 
emplified by the Omori law of aftershocks and the Gutenberg- 
Richter relation for the number of earthquakes above a given 
magnitude, and more recently, by the fractal properties of 
earthquake spatial occurrence ll|,|2|,|3j,|4|]. Less attention how- 
ever has been paid to the timing of individual earthquakes, for 
which a unifying picture was missing until the work of Bak et 
al. |5J,|g,|7|]. The main relevance of that work was the sem- 
inal introduction of scaling concepts in earthquake statistics, 
which provide a powerful tool to unify descriptions and to de- 
rive relations between different quantities (in this case, as we 
will see, by scaling we do not mean just power-law relations, 
as it is sometimes assumed). Later, Bak el a/.'s procedure was 
modified to study the distribution of times between consecu- 
tive events in a single spatial region |8|. 

Let us consider the seismicity of an arbitrary spatial region. 
Given a lower bound M c for the magnitude, and intentionally 
disregarding the spatial degrees of freedom, a marked point 
process in time of the form (to, Mo), (t\,M\), ... is obtained, 
where f, denotes the time of occurrence of event i, with a mag- 
nitude Mi > M c . The recurrence times are defined as the time 
intervals between nearest-neighbor (i.e., consecutive) events, 
%i = tj — ti- 1 . In the case of stationary seismicity (character- 
ized by a long-term linear relation between the accumulated 
number of earthquakes and time), for spatial regions of linear 
size ranging from 20 km to the whole world, and for magni- 
tude bounds from 1.5 to 7.5, the probability densities D{x) of 
the recurrence time were found to verify a universal scaling 

law Hi, 

D(x)=Rf(Rx), (1) 

where / is a universal scaling function and the scaling factor R 
is the rate of seismic occurrence, defined as the mean number 
per unit time of events with M > M c in the region, and given 
by the Gutenberg-Richter law, R(M C ) = R \Q~ bMc , with the 
b— value usually close to 1 . 

As no separation of mainshocks and aftershocks is per- 
formed, the recurrence-time distribution consists of a mixture 
of different aftershock sequences and more or less indepen- 
dent events; therefore, it is very surprising that distinct regions 
and earthquakes of disparate sizes present such a extreme de- 
gree of regularity. 



But even more surprising, the scaling relation for the 
recurrence-time distribution reveals that seismicity is in a 
highly orchestrated state, in which the removal of events 
(when the lower bound M c is raised) does not affect the prop- 
erties of seismic occurrence, as the distribution keeps the same 
shape (with only a different mean) independently of M c . In 
general, when some events are removed from a point process, 
the properties of the process do change; therefore, the dis- 
tribution of recurrence times constitutes a very special case, 
invariant under a transformation akin to those of the renor- 
malization group (RG) in real space |3, 9, 10]. 

The first step of our renormalization-group transformation 
consists on the raising of the lower bound M c . This im- 
plies that only a fraction of events survives the transformation, 
which defines a different recurrence-time distribution. (If M c 
is increased in one unit, we are dealing with an authentic dec- 
imation, as only about 1 tenth of events are kept, due to the 
Gutenberg-Richter law.) The second part of the procedure 
is the scale transformation, which changes the time scale to 
make the new system comparable with the original one. 

The Poisson process, characterized by an exponential 
recurrence-time distribution, represents a trivial solution to 
this problem when there are no correlations in the process and 
therefore events are randomly removed. Indeed, it has been 
argued that the scaling function / can only be an exponential 
function II ill . However, the scaling function clearly departs 
from an exponential |8], and in this way the relevance of cor- 
relations in the structure of seismicity becomes apparent. 

Indeed, the scaling function / can be described by a de- 
creasing power law for intermediate times, Rx < 1, with an 
exponent about 0.3, and a faster decay for long times, Rx > 1, 
well approximated in both cases by a gamma distribution. No 
model of earthquake occurrence is assumed to obtain these re- 
sults, they are a fundamental characteristic of seismicity. For 
short times, Rx < 0.01, the condition of stationarity is usually 
not fulfilled and the behavior is not universal. Nevertheless, in 
the non-stationary case the process can be transformed into a 
stationary one with an appropriate transformation of the time 
axis, and then the same scaling relation is found to hold again 

It should be noted that the term "correlations" can be un- 
derstood in two forms: If the recurrence-time distribution is 



2 



not an exponential, this implies the existence of a memory ef- 
fect in the process, as events do not occur independently at 
any time, as it would be the case for a Poisson process. But 
further, the recurrence time may depend on the history of the 
process, in particular the occurrence time and magnitude of 
previous events. We shall see that this type of correlations 
are responsible of the breaking of the memoryless character 
of seismicity. 

In general, the time between two consecutive earthquakes, 
T;, may depend on the magnitude of the previous event, 
Mp re = Mi—\, the previous recurrence time, T/_i, and also on 
the occurrence of preceding events, i — 2, i— 3, etc. In then- 
turn, the magnitude of the i— th event M, can depend on t,, 
Aff_ i , t,-_ i , and so on 1 1 2 ] . We shall only consider here the de- 
pendence of Xi with M,_ i , as we have verified it is the most im- 
portant (together perhaps with the dependence of T, with T,-_ i ). 
Note also that although the dependence of the recurrence time 
and the magnitude with the distance between events can be 
important, as we have not considered spatial degrees of free- 
dom, we do not need to take this effect into account. 

So, in what follows we study the effect in the structure 
of seismicity of the simplified case of correlations between 
the recurrence time and the magnitude of the previous earth- 
quake. If we raise the magnitude threshold from M c to M' c the 
distribution of recurrence times for events with magnitudes 
M > M' c can be obtained from the distribution for events with 
M > M c . Assuming that an event with magnitude Mq > M' c 
has occurred, we can write for the next event above (or at) M' c , 

P[ recurrence time > x for events M > M' c ] = 
P[ 1st return time > x and Mi > M' c ] + 
P[ 2nd return time > x and M\ < M' c and M 2 > M' c ] + 
P[ 3rd return time > T and Mi <M' C ,M 2 < M' c , and M 3 > M' c ] 

(2) 

where the sum has to be extended up to infinity. P denotes 
probability and the n— th return time is defined, for events with 
M > M c , as f, — ti- n , that is, as the elapsed time between any 
event and the n— th event after it (of course, the 1st return time 
is the recurrence time). As the recurrence time depends only 
on the previous magnitude, but the magnitude is independent 
on any other variable we can write 

P[ recurrence time > x for events M > M' c ] = 
P[ 1st return time > x \M\ > M[]p+ 
P[ 2nd return time > x \M\ < M' C ,M2 > M' L ]qp+ 
P[ 3rd return time > x\M\ < M' C ,M 2 < M' cl Mi, > M' c ]q 2 p - 



with 



p = P[M >M' C \M> M c ] = IQ-Wc-Mc) j 



(3) 



(4) 



where the last equality comes from the Gutenberg-Richter law 
and q = 1 - p = P[M <M' C \M> M c }. 

Derivation in Eq. (0 with respect x yields the probabil- 
ity densities D of the different return times; as the recurrence 
times are independent on each other, we use that the n— th- 
return-time distribution is given by n convolutions of the first- 



return-time distributions (denoted by the symbol *) to get 

T 1/2 D(T) = P D 1 (z)+qpD^x)*D l (x)+q 2 pD 1 (x) *D±(x) *D { (t) - 

(5) 

where ~T ii 2 D{x) denotes the probability density for events 
with M > M' c as a transformation T^n of the probability den- 
sity for events with M > M c , D(x); more precisely, D(x) = 
D(x\M>M c ), andT 1/2 D(T) =D(x\M>M' c ). The subscript 
1 /2 refers to the fact that this is only the first half of the RG 
transformation. D| and D i denote the recurrence-time prob- 
ability densities for events above M c conditioned to the fact 
that the magnitude of the previous event is above or below M' c 
(t or I), respectively. To be precise, D^(x) = D(x\M pre > 
M' C ,M > M c ), and Dj_(t) = D(x \M pre < M' C1 M > M c ). 
In Laplace space the things are simpler, 

T 1/2 D(s) - pDj (s) + qpDj (s)D l (s) +q 2 P D ] (s)D l (s)D l (s) + ■■■ 

(6) 

Notice that we have used the same symbol D for both the 
probability densities and for their Laplace transforms (which 
we may call generating functions), although they are differ- 
ent functions, of course. As q and £>|(s) are smaller than one 
(this is general for generating functions), the infinite sum can 
be performed, turning out that 

p£> T (s) pD^s) 



1-qD^s) l-D(s)+pD^ S y 



(7) 



using that D(x) is in fact a mixture of the distributions Dj and 
D;, of the form D = pDj + qD±. 

The previous equation Q describes the first part of the 
transformation. The second part is the scale transformation, 
which puts the distributions corresponding to M c and M' c on 
the same scale. We will obtain this by removing the effect of 
' the decreasing of the rate, which has to be proportional to p, 
so, 

T l/2 D(x)^p- x T x/2 D{x/p) : (8) 

and in Laplace space we get 

T 1/2 D(j) — > Ti /2 D(ps). (9) 

Therefore, the combined effect of both transformations leads 
to 

pD^ps) 



TD(s) 



l-D(p S )+pD 1 {ps)- 



(10) 



In order to get some understanding of this transformation 
we can consider first the case in which there are no correla- 
tions between the magnitude and the subsequent recurrence 
time. Then, Dt = = D = Do and 



TDo(s) 



pD (ps) 



(11) 



i - qD (ps) ' 

The fixed points of the transformation are obtained by the so 
lutions of TDq(s) = Dq(s). If we introduce © = ps and sub 
stitute p = co/s, q = 1 — o/s, we get, separating variables, 



1 



1 



1 



sDo(s) 



(oD ((o) 



1 

— =k; 

CO 



(12) 
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where A; is a constant, due to the fact that p and s are indepen- 
dent variables and so are s and co. The solution is then 



Do(s) 



1 



l + ks' 



(13) 



which corresponds to the Laplace transform of an exponential 
distribution. Indeed, 



A>(t) 



-t/k 



(14) 



So, in the case in which there are no correlations in the pro- 
cess, the only scale invariant distribution is, as one could have 
expected, the exponential distribution, characteristic of the 
Poisson process. Let us see how the existence of correlations 
between the magnitudes and the recurrence times changes this 
picture. 

Correlations introduce new functions in the process. In our 
case, in order to iterate the transformation T we need to know 
how D-f transforms as well. It turns out that Dj verifies an 
equation very similar to Eq. (1101 . which depends on Df (t) = 
D(x \M P re > M",M > M c ). So, in order to apply again T we 
also need an equation to transform Da, which in turn depends 
on higher values of the magnitude threshold. In this way we 
obtain a hierarchy of equations. An easy way to break this 
hierarchy is to assume that, at least at the fixed point, Df has 
the same functional form as D, but in a different scale. So, let 
us assume 



D t (t)=AD(At), 



(15) 



where A depends on M' c -M c with A(0) = 1. 

Figure [2] illustrates this behavior using worldwide earth- 
quakes from the NEIC PDE catalog |8]. The distributions Z?j 
for different values of M' c keeping M c = 6 collapse onto a sin- 
gle curve under rescaling of the axes. For each distribution the 
scaling factor is the inverse of its mean value, /?t, and there- 
fore A = R^/R. The behavior of A as a function of M' c appears 
in Fig. |2 A flat line would indicate absence of correlations, as 
Rt would be identical to R and therefore Dt = D. In real seis- 



micity, R^ increases with the magnitude of the previous event, 
which means that the mean time between events decreases. In 
the figure, fits of the type A(M' C -M c ) =A+ C{M' C - M c ) and 
A(M' C -M c ) = Ae c( - M '^ M ^ are shown; in both cases it turns 
out to be that A ~ 1 and C is in the range 0.18-0.20. 

Returning to our calculation, in Laplace space D^(s) = 
D(s/A); therefore, the transformation (II 01 turns out to be 



TD(s) 



pD(ps/A) 



l-D(ps)+pD(ps/A)' 



(16) 



As this discrete transformation is difficult to deal with, we 
will look at the infinitesimal transformation defined by M' c — > 
M c . Introducing 8 =M' C ~M C , this implies p = \Q- b ( M 'c- M c) ~ 
1 B8 with B = MnlO, D(ps) = D(s) - BsD'{s)8, A ~ 1 + 
C<5, and D(ps/A) = D(s) -(B + C)sD'(s)8. Substituting in 
Eq. dl6l and up to first order in 8 we get 



TD(s)~D(s)+{M*)- 



l][BD(s)+CsD'(s)]-BsD'(s)}8. 

(17) 



In order to find the fixed point of the infinitesimal transforma- 
tion we impose that the coefficient of 8 is zero, obtaining, 



D'(s) 



BD(s)[l-D(s)] 



s{B + C[l-D(s)}Y 
whose integration yields to 

ksD 1+£ {s)+D(s)- 1 =0, 



(18) 



(19) 



where the exponent 1 + e comes from the definition e = C/B. 

We immediately see that in the case of no correlations, 
C = 0, and therefore e = 0, recovering Eq. Jl 3i and then the 
exponential form for D(x). But there are other values of e for 
which the previous equation can be easily solved. 

Let us consider first the case 8 = 1, corresponding to B = C. 
The solution of Eq. dl9l is 



D{s) 



^f\TWs- 1 
2ks ' 



(20) 



whose inverse Laplace transform can be calculated [13], turn- 
ing out to be, 



D(T) 




(21) 



with erf the error function. The asymptotic behavior of 
D(x) is very clear: for x — > it diverges as a power law, 
D(x) — > 1/ V nkx, whereas for x — > °°, D(x) decays exponen- 
tially, using the expansion of the error function 1 13] and the 
fact that a power-law factor varies much more slowly than an 
exponential, for large arguments. 

It is interesting also to study the case e — > 0, characteristic 
of weak correlations, C ~ 0. We can write the solution D(s) as 
a perturbation of the Poisson behavior corresponding to e = 0, 
i.e., D(s) = Dq(s) +u(s), with D (s) = (1 +ks)~\ see Eq. 
( fT3l . Substituting into Eq. ([19} and using (1 +ks) £ ~ 1 + 
eln(l +ks) we get 



(1+ksY 



1 



1 +ks 



(22) 



and therefore D(s) is just a power of the transform of the ex- 
ponential density, which means that D(x) is a gamma distri- 
bution, i.e., 



D(x) 



1 



fcT(l-e) \x 



-T/k 



(23) 



In this way, for e ~ 0, we also get a power law for short times 
and an exponential decay for long times. 

This behavior is by no means exclusive of e = 1 or e ~ 0. 
If in Eq. ( 1191 we consider the limit s^<»we get, as D(s) — * 
(which is general for generating functions), 



D(*)->l/(fa)TO, 



(24) 
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and, if 1 + e > 0, by means of a Tauberian theorem 1 14], 



D(T) 



1 



*r[i/(i+c)] \t 



, for T — > 0. 



(25) 



This implies that for all e > 0, which corresponds to A > 1, 
and then to a shortening of the recurrence times after large 
earthquakes, we get a decreasing power law, which is a signa- 
ture of clustering, and in concordance with real data |8]. (The 
situation would be reversed if —1 < e < 0.) 

Further insight can be obtained from Eq. dl9> . writing it 
as D = 1 — ksD l+£ , which can be iterated to get the Taylor 
expansion of D(s). Up to second order, it turns out that D(s) = 
1 — ks + (1 + e)k 2 s 2 , and as the coefficients of s n yield the 
moments of the distribution (with a (— l)"/n! factor (3) we 
conclude that (t) = k, (t 2 ) = 2 (1 + e)k 2 and the coefficient of 
variation is cv = \J (t 2 ) - (t) 2 / (t) = Vl +2e. 

In fact, the Taylor expansion can be obtained up to any or- 
der, and so, all the moments (t") exist and D[x) decays faster 
than any power law for T — > °°, in agreement again with the 
observations [8]. 

To conclude, we note that both B and C can be understood 
as critical exponents. As the energy radiated by an earthquake 
is about E = £bl0 3M / 2 0]> the Gutenberg-Richter law writes 
P[E > E c ] °c l/Ec b/3 . Also, assuming the exponential form 
for A, frt = AR = R e c(K-M c ) and we get 



'T 



2C 
JB i() 



(26) 



If we define v = 2C/(31nlO), this exponent, the b— value, 
and the exponent of the recurrence-time distribution for short 
times given by Eq. \25\ fulfill the following scaling relation, 

e C v 



1+6 B+C 



-2^/3' 



(27) 



just using £ = C/B and B = b In 10. In fact, we must under- 
stand this relation as the contribution of the M pre — t correla- 
tions to the recurrence-time distribution, as the value obtained 
for C (about 0.2) is too small to account for the value of the 
exponent, about 0.3, using a b— value ~ 1. We believe the in- 
clusion of other correlations in the calculation will yield to a 
better quantitative agreement. 

Our approach mainly consists of a simplification of real 
seismicity, which allows to understand the complex structure 
of seismic occurrence in the time and magnitude domains. It 
is remarkable that simply by imposing the self-similarity of 
the process and with the only assumption of the scaling of the 
conditional distributions Dj (which is in agreement with the 
observations), we get such a general characterization of the 
recurrence-time distribution. With this study we have shown 
how the structure of seismic occurrence in time, space, and 
magnitude can be understood as a critical phenomenon and 
then constitutes a statistical -physics problem 1 15]. 
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FIG. 1 : Recurrence- time probability densities conditioned to M pre > M' c , M' c varying from 6 to 7 and with M > 6, for worldwide earthquakes 
from 1973 to 2002. The data collapse indicates that Dt (t) verifies the same scaling relation as D(x). The straight line is a fit to D(t), turning 
out to beoc e -Vi-4/T M . 
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FIG. 2: Inverse of the mean recurrence time, R+, scaled by R for recurrence periods started by events with M pre > M' c and ending with M > 5. 
The data correspond to worldwide earthquakes from 1973 to 2002. The error bars mark two standard deviations of the mean value, and the 
two curves are the linear and exponential fits explained in the text. The last two points have not been taken into account for the fits. 



